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METHOD FOR PROVIDING ROBUSTNESS TO THE SOLUTION OF THE 

CUBIC EQUATION OF STATE 

1 Field of the Invention 

This invention relates to models which are used in 
process simulation and more particularly to the cubic 

equations of state. 

2 Description of the Prior Art 

The use of computer simulations to develop new 
processes, design and retrofit plants, and optimize the 
operation of existing plants is widely practiced in the 
process industry. Throughout process design, 

development, and operation, the mathematical models that 
are used in process simulation can be used to make better 
engineering and business decisions. Rigorous models 
require complex numerical solvers. The stability of the 
numerical solvers can be improved by using integrated 
physical property calculation solvers. 

The cubic equations of state (CEOS) relate pressure, 
temperature and molar volume of a mixture. The deviation 
of the physical properties of a mixture from ideality is 
expressed by the compressibility factor. The 
compressibility factor is calculated from the CEOS, where 
the compressibility factor is a function of volume terms 
raised to third power. The CEOS offers the advantage of 
using very few standard physical properties (critical 
pressure, critical temperature, critical compressibility 
factor, critical volume, and accentric factor) as input 
to evaluate the actual compressibility factor. 

All other physical property equations are functions 
) of the compressibility factor. These include fugacity 
enthalpy, entropy, viscosity, volume, density and 
internal energy. The most commonly used forms of CEOS 
are Soave-Redlich-Kwong (SRK) , Peng-Robinson (PR), 
Redlich-Kwang (RK) , and Ideal equations of state. 

The solution to the CEOS to compute the 
compressibility factor encounters three problems: 
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non-unique roots; 

• trivial roots; and 

• singularities. 

The two phase Vapor-Liquid equilibrium region is 
discontinuous in nature. The solution of the CEOS in a 
two phase Vapor-Liquid equilibrium region can sometimes 
result in three non-unique roots. In a narrow region of 
physical conditions, the CEOS can give one valid root for 
the liquid phase and one valid root for the vapor phase. 
Under most conditions of pressure, temperature, and 
composition, however, no root can be found either for the 
vapor phase or the liquid phase. This is known as a 
trivial roots problem and is commonly encountered over a 
wide range of operating conditions. The third numerical 
problem, singularity, arises when the valid root happens 
to be at a local maximum or minimum of a pressure-volume 
curve. This represents infinite possible values of 
derivatives of dependent properties such as K value, 
enthalpy departure, and fugacity coefficients. 

Physical property prediction tools have overcome 
these deficiencies by identifying vapor, liquid and 
vapor-liquid regions to determine the conditions 
necessary to identify solution of the CEOS. However, 
these procedures can sometimes not be particularly 
robust, especially when a system either approaches or 
experiences a discontinuity during the course of the 
solution . 

The prior art property prediction tools are 
described in: 

1. Poling B. E. , Prausnitz J. M., O'Connell J. P., The 
Properties of Gases & Liquids, 5 th edition, McGraw- 
Hill, November 2000 ("Poling et al."). 

2. Agarwal R, Li Y. K. , Santollani 0., Satyro M. A., 
Vieler A., Uncovering the Realities of Simulation 
Part II, Chemical Engineering Progress, pp 64-72, 
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June 2001 ("Agarwal et al."). 
3 Twu C. H., Sim W. D., Tassone V., Getting a Handle on 
Advanced Cubic Equations of State, Chemical 
Engineering Progress, pp. 58-65, November 2002 ("Twu 

The present invention defines a procedure whereby 
the solution of the CEOS can be represented in continuous 
form overcoming many of the robustness issues experienced 
by the prior art. In arriving at the solution of the 
10 present invention it was determined that the problems in 
the solution of the CEOS, namely, non-unique roots, 
trivial roots and singularities were primarily 
responsible for the failure of the solver. Thus the 
method of the present invention was developed to 
eliminate all of these problems and generate a continuous 
solution for compressibility factors over the entire 
range of operating conditions. 

Summar y, of the Invention 
A method for calculating the compressibility factors 
20 of a mixture using the cubic equation of state (CEOS). 
The mixture has a liquid region. The method comprises: 

a) calculating a pure liquid compressibility 
factor when the compressibility factor is in the liquid 

region; _ . c 

b) using the second virial coefficient of the CEOS 

to determine from the pure liquid compressibility factor 
a vapor compressibility factor extensible in the liquid 

region; and _ 

c) accepting the higher of the pure liquid 
30 compressibility factor or the vapor compressibility 

factor extensible in the liquid region. 

A method for calculating the compressibility factors 
of a mixture using the cubic equation of state (CEOS). 
The mixture has a vapor-liquid equilibrium region. The 

method comprises: 

a) determining that the compressibility factor is 
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in the vapor-liquid equilibrium region; 

b) calculating when the compressibility factor is 
in the vapor-liquid equilibrium region a pure liquid 

compressibility factor; 

c) using the virial coefficients of the CEOS and 
the pure liquid compressibility factor calculated when 
the compressibility factor is in the vapor-liquid 
equilibrium region to determine a vapor compressibility 
factor in the vapor-liquid equilibrium region; and 

d) accepting the higher of the pure liquid 
compressibility factor calculated when the 
compressibility factor is in the vapor-liquid equilibrium 
region or the vapor compressibility factor calculated in 
the vapor-liquid equilibrium region. 

A method for calculating the compressibility factors 
of a mixture using the cubic equation of state (CEOS). 
The mixture has a vapor-liquid equilibrium region. The 

method comprises: 

a) calculating a pure liquid compressibility 

factor for the vapor-liquid equilibrium region; 

b) using the virial coefficients of the CEOS and 
the pure liquid compressibility factor calculated for the 
vapor-liquid equilibrium region to determine a vapor 
compressibility factor in the vapor-liquid equilibrium 
region; and 

c) accepting the higher of the pure liquid 
compressibility factor calculated for the vapor-liquid 
equilibrium region or the vapor compressibility factor 
calculated in the vapor-liquid equilibrium region. 

A method for calculating the compressibility factors 
of a mixture using the cubic equation of state (CEOS). 
The mixture has a liquid region and a vapor-liquid 
equilibrium region. The method comprises: 

determining that the compressibility factor is 
either in the liquid region or in the vapor-liquid 
equilibrium region; 
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the method when the compressibility factor is in the 

liquid region comprising: 

calculating a pure liquid compressibility factor; 

using the second virial coefficient of the CEOS to 
determine from the pure liquid compressibility factor a 
vapor compressibility factor extensible in the liquid 
region; and 

accepting the higher of the pure liquid 
compressibility factor or the vapor compressibility 
factor extensible in the liquid region; 

the method when the compressibility factor is in the 
liquid region comprising: 

calculating a pure liquid compressibility factor; 

using the virial coefficients of the CEOS and the 
pure liquid compressibility factor to determine a vapor 
compressibility factor in the vapor-liquid equilibrium 
region; and 

accepting the higher of the pure liquid 
compressibility factor or the vapor compressibility 
factor calculated in the vapor-liquid equilibrium region. 

A method for calculating the compressibility factors 
of a mixture using the cubic equation of state (CEOS) . 
The mixture has a liquid region and a vapor-liquid 
equilibrium region. The method comprises: 

a) calculating a pure liquid compressibility 
factor when the compressibility factor is in the liquid 
region; 

b) using the second virial coefficient of the CEOS 
to determine from the pure liquid compressibility factor 
a vapor compressibility factor extensible in the liquid 
region; and 

c) accepting the higher of the pure liquid 
compressibility factor or the vapor compressibility 
factor extensible in the liquid region; 

d) calculating a pure liquid compressibility 
factor for the vapor-liquid equilibrium region; 
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e) using the virial coefficients of the CEOS and 
the pure liquid compressibility factor calculated for the 
vapor-liquid equilibrium region to determine a vapor 
compressibility factor in the vapor-liquid equilibrium 
region; and 

f) accepting the higher of the pure liquid 
compressibility factor calculated for the vapor-liquid 
equilibrium region or the vapor compressibility factor 
calculated in the vapor-liquid equilibrium region. 

A method for calculating the liquid compressibility 
factor of a mixture using the cubic equation of state 
(CEOS). The mixture has a liquid region. The method 
comprises : 

a) calculating a pure liquid compressibility 
factor when the compressibility factor is in the liquid 
region; 

b) determining that the rate of change of the pure 
liquid compressibility factor with temperature is or is 

not less than zero; and 

c) accepting the pure liquid compressibility as 
the liquid compressibility factor when the pure liquid 
compressibility factor rate of change is not less than 
zero . 

A computer readable medium having instructions for 
causing a computer to execute a method for calculating 
the compressibility factors of a mixture using the cubic 
equation of state (CEOS), the mixture having a liquid 
region. The instructions comprise: 

a) calculating a pure liquid compressibility 
factor when the compressibility factor is in the liquid 
region; 

b) using the second virial coefficient of the CEOS 
to determine from the pure liquid compressibility factor 
a vapor compressibility factor extensible in the liquid 
region; and 

c) accepting the higher of the pure liquid 
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compressibility factor or the vapor compressibility 
factor extensible in the liquid region. 

A computer readable medium having instructions for 
causing a computer to execute a method for calculating 
the compressibility factors of a mixture using the cubic 
equation of state (CEOS), the mixture having a vapor- 
liquid equilibrium region. The instructions comprise: 

a) calculating a pure liquid compressibility . 
factor for the vapor-liquid equilibrium region; 

b) using the virial coefficients of the CEOS and 
the pure liquid compressibility factor calculated for the 
vapor-liquid equilibrium region to determine a vapor 
compressibility factor in the vapor-liquid equilibrium 
region; and 

c) accepting the higher of the pure liquid 
compressibility factor calculated for the vapor-liquid 
equilibrium region or the vapor compressibility factor 
calculated in the vapor-liquid equilibrium region. 

A computer readable medium having instructions for 
causing a computer to execute a method for calculating 
the compressibility factors of a mixture using the cubic 
equation of state (CEOS), the mixture having a liquid 
region. The instructions comprise: 

a) calculating a pure liquid compressibility 
factor when the compressibility factor is in the liquid 
region; 

b) determining that the rate of change of the pure 
liquid compressibility factor with temperature is or is 

not less than zero; and 

c) accepting the pure liquid compressibility as 
the liquid compressibility factor when the pure liquid 
compressibility factor rate of change is not less than 
zero . 

Description of t he Drawing 
Fig. 1 shows a typical pressure (P) - volume (V) 
relationship for a hydrocarbon compound. 
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Fig 2 is the flowchart for the method of the 
present Invention to calculate the vapor compressibrlrty 

^^"Fig. 3 is the flowchart for the method to calculate 
the liquid compressibility factor. 

Fig 4 shows in block diagram form a process plant 
in which the method of the present invention can be used. 
DescriEtion " f the Pre f erred Embodrment (s) 
Referring now to Fig. 1 there is shown by a dotted 
line a typical pressure <P, - volume (V) relationship for 
a hydrocarbon compound. The P-V relationshrp has, as is 
described in Poling et al., three distinct * 
high pressure region 10 representing the Irqurd phase, 
shows a sharp variation. A low pressure 12 regron 
representing the vapor phase shows a slow varratron An 
intermediate region 14 representing vapor-Uqurd 
equilibrium region shows a moderate varratron wrth a 
complex behavior. The P-V relationship in the ,^or-li^ 
equilibrium region is not smooth and can have local maxrma 
14a and minima 14b as is shown in Frg. 1. 

While the techniques of the prior art descrrbed rn 
Poling et al., Agarwal et al. and Twu et al. can as rs 
shown by the dashed lines in Fig. 1 closely aPP™te 
the pure vapor and pure liquid regions they devrate from 
the actual data in the vapor-liquid equilibrium regron. 
The prior art has tried several approaches to overcome 
this deviation, including extension of pure vapor rn to 
iiguid phase and extension of pure liquid in to vapor 
phLe. These approaches fail to approximate the vapor 
0 liquid region as the prior art solutions are almost always 
unstable and inaccurate in that region. 

The methods of the present invention rmplemented for 
solving the CEOS for the vapor and liquid compressibilrty 
factors can as is shown in Fig. 1 by the solrd Irne 
closely approximate the pure liquid and vapor regron as 
well as the vapor-liquid equilibrrum regron. These 
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methods are described below. The method of the present 
invention is focused on the P-V prediction in the vapor- 
liquid equilibrium region by : 

(a) improving stability 

(b) improving accuracy 

Fig. 2 is the flowchart for the method to calculate 
the vapor compressibility factor and Fig. 3 is the 
flowchart for the method to calculate the liquid 

compressibility factor. 

As is described in Poling et al. and Agarwal et al., 
the CEOS that is solved to obtain the value of Z, the 
compressibility factor, is 

(E. 1) 

Z 3 +B*Z 2 +C*Z + D = 0 

where the parameters B, C and D are known as second, 
third and fourth virial coefficients. For example, under 
the Soave-Redlich-Kwong system, their definitions are as 

follows: 

B = -l 

C = A'-B'-B' 2 
D = -A'* B' 

Where 

aP 



A = 



B' = 



R 2j2 

bP 



30 



RT 

CALCULATIONS FOR Z r , VAPOR COMPRESSIBILITY FACTOR 

Once the B, C, and D values are computed, the solver 
goes through the sequence of equations and conditions 
shown in Fig. 2 and described below to determine the 
compressibility factor. The sequence of equations and 
conditions followed by the solver depends on the value of 
R The method follows the CEOS solution in the vapor 
region. It uses liquid compressibility as an extensron 
of vapor compressibility in the liquid region. The 
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equations described below are used to smooth out the Z v 
calculations during transition from vapor to vapor-liquid 
to liquid region. 

Equations Used In Calculation Of Z v 
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where the value of the variable R determines the possible 

solutions of the CEOS. 

When R > 0, the CEOS has only one real root and two 
conjugate complex roots. The solution is either in a 
pure vapor or pure liquid region. In the vapor region, 
the compressibility factor for pure vapor is accepted. In 
the liquid region, the compressibility factor for pure 
liquid is accepted as an extension of the vapor 
compressibility factor. It should be noted that it is a 
common practice to hypothetically extend the pure vapor 
solution in to the pure liquid region. However, this can 
cause a sharp change in compressibility factor leading to 

instability in the solver. 

When R ^ 0, the CEOS can have either two or three 
roots. When R = 0 there are three real roots and at least 
two of them are equal. When R < 0 the solution is three 
unequal roots. The solution falls in the vapor-liquid 
equilibrium region. This is where state changes occur 
and the solution of the CEOS becomes discontinuous The 
solution faces, as was described above, problems of non- 
unique roots, trivial roots and singularities. 

This overall scheme is elaborated in following 

discussion. 

r- (21*D-9*B*C + 2*B 3 } (V.2a) 
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(V.3a) 



where Z L is the liquid compressibility factor. 



7 =3 
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r 21*D-9*B*C + 2*B 
I 54" 



3\ 



21*D-9*B*C+2*B 



3^ 



54 



(V. 3b) 



5 + Z t 



^g+Z >#+^); -4*(C+Z £ *(B + Z L )) 

2 



(V.4) 



(V.5) 



(V.6) 



D 



„ here in equation (V.6). fc i- the pure liquid 
compressibility factor for the vapor-liquid equrlrbrrum 
region determined by equation (V2.b). 

Mathematics, see for example Editors : Perry R. H., 
Green D . Maloney J. 0., Perry's Chemical Engineers 

Handbook, 7* edition, McGraw-Hill, July 1999 has 
established that when R * 0, equation (E.l) has three 
roots 



Z = P + Q 



(E.2) 
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r+Q . - Q ) (E ' 3) 



Z = r^- + 



P + Q _iJ[( p _n) 



(E.4) 



2 2 



Where 



i 2 = -1 



P = Jzl{21*D-9*B*C + 2*B 1 }+jR 



(E.5) 



g = 3 |zi{27*D-9*B*C + 2*5 3 }-V^ (E.6) 

As is shown in Fig. 2, the method of the present 
invention starts with the calculation of R using equation 
(V 1) If R > 0, there is, as is shown in equations E.2 
to' E 4, only one real root and two complex conjugate 
roots. The one real root represents a solution either in 
the pure vapor region or in the pure liquid region. 

The method proceeds to equation (V.2a) where it 
determines if the square root of R is greater than a 
certain mathematical value. If this is false then the 
compressibility factor is in a vapor region. The vapor 
compressibility factor value is accepted as the pure 
vapor compressibility factor determined from equation 

(V.3b) . , . 

If the evaluation of equation (V.2a) determines that 

the square root of R is greater than a certain value then 

the compressibility factor is in a liquid region. The 

method proceeds first to equation (V.3a) to obtain the 

pure liquid compressibility factor. It then applies 

correction as in equation (V.4) to obtain the vapor 

compressibility factor and finally it accepts the higher 
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of the two values, V.3a or V.4, as actual vapor 
compressibility factor as in equation (V.5). 

If R = 0 then all three roots of CEOS are real. Out 
of these three roots, at least two are equal. This 
condition is however lumped with the condition R < 0 to 
avoid discontinuities. 

When R < 0, then all three roots given by equations 
(E.2), (E.3), and (E.4) are real and unequal. This 
represents a vapor-liquid equilibrium region. But the 
roots given by equations (E.2), (E.3), and (E.4) are 
unrealistic. Instead, the roots given by equation (V.2b) 
are mathematically accepted. Corresponding liquid and 
vapor compressibility factors are evaluated by equations 
(V.2b) and (V.6) respectively. The greater of these two 
values, V.2 or V.6, is accepted as a true vapor 
compressibility factor for the vapor-liquid equilibrium 

region . 

CALCULATIONS FOR Z L , LIQUID COMPRESSIBILITY FACTOR 

The algorithm follows the CEOS solution in the liquid 
region. It uses vapor compressibility as the extension of 
liquid compressibility in the vapor region. Equation 
V.2b is used to smooth out the Z t calculations during the 
transition from vapor to vapor-liquid to liquid region. 
Again the logic behind this framework of equations, the 
flowchart for which is shown in Fig. 3, is same as that 
discussed for vapor compressibility factor calculations. 



30 



dZ 



2 dB ^ dC dD 
Z , 2 — + Z L — + 



dT 



dT dT 
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Ll ~ dC (L.3) 
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When R < 0 in the vapor-liquid equilibrium region, the 
compressibility factor given by equation (V.2b) is used. 

When R > 0 the solution lies either in the pure 
liquid region or the pure vapor region. Equation V.2a is 
determined and if it is true then the pure liquid phase 
is assumed. Otherwise if equation V.2b is false then the 
pure vapor phase is assumed. 

In the pure liquid phase, first the liquid 
compressibility factor is determined using equation 
(V 3a) . Then the local analytical derivative of liquid 
compressibility factor is determined using equation 
(L.l). If the derivative is * 0 then the pure liquid 
compressibility factor is accepted as it is. If the 
derivative is less than zero then a modified liquid 
compressibility factor given by simplified equation (L.3) 
is used to avoid the unrealistic solution region. 

in the pure vapor phase, first the vapor 
compressibility factor given by equation (V.3b) is 
determined. Equation (L.2) then uses this vapor 

compressibility factor to evaluate a pseudo liquid 
compressibility factor as an extension of the vapor 

compressibility factor. 

The Soave-Redlich-Kwong, Peng-Robinson, Redlich- 

Kwang, and Ideal forms of the CEOS have all been used to 

validate the method of the present invention for solving 
the CEOS for the compressibility factor. 

It should be appreciated that the method of the 
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present invention can accelerate the calculations of 
physical property by using ideal property values as 
initial guesses for iterative calculations, reusing 
earlier calculation results, either intermediate or 
final, after they have been generated. Further two types 
of caches to hold previous calculation results, namely, 
the intermediate results cache and the final results 
cache can be designed. With the intermediate results 
cache, once a quantity such as a compressibility factor 
has been generated to calculate a certain physical 
property, it can be reused for evaluating other 
properties for the same inputs. With the final results 
cache, a property has been calculated through linear 
extrapolation using the value of the property and its 
derivatives that have been obtained before for delta 

change in conditions. 

Appropriate scaling factors are incorporated to 
account for different variations in mathematically 
related values. For example, the exponential 

relationship between vapor pressure and temperature is 
accounted for, by using natural log of vapor pressure 
while linearly changing temperature during iterative 
procedures . 

As can be appreciated by those of ordinary skill in 
the art, the present invention may be used in connection 
with any industrial process that uses gases, liquids and 
their mixtures. Examples of such processes include, 
without limitation, those related to oil & gas, chemical, 
petroleum and pulp & paper. As is shown in Fig. 4, the 
industrial process 20 usually resides in a plant 22 which 
also has an automation system 24 for measurement and 
control of the process 20 in a manner well known to those 
in the art. The method of the present invention may 
reside on a computer readable medium and thus can be 
loaded into a computing device 26 such as a personal 
computer. The computing device 26 may, as is shown by 
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the dotted lines, be associated with the industrial 
process 20 if only measurements are desired or a part of 
or associated with the automation system 24 if both 
measurement and control is desired or associated with 
both process 20 and system 24. 

It is to be understood that the description of the 
preferred embodiment (s ) is (are) intended to be only 
illustrative, rather than exhaustive, of the present 
invention. Those of ordinary skill will be able to make 
certain additions, deletions, and/or modifications to the 
embodiment ( s ) of the disclosed subject, matter without 
departing from the spirit of the invention or its scope, 
as defined by the appended claims. 



